# Load necessary libraries
library(ggplot2)
library(dplyr)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
library(viridis)
library(doBy)
library(gganimate)
library(ggplot2)
library(maps)
library(ggthemes)

#Set wd for the plots
setwd("~/OneDrive - Indiana University/FromGoogle/syd/plots/")
# Collapse data by gid (mean of variables) and calculate the outbreak rate
collapsed_data.1 <- summaryBy(clear_zoon_confirm+NDVI_mean+fires.cell.count~latitude+longitude,
                              data=full.2003.2018, FUN=c("mean"), na.rm=T)
collapsed_data.2 <- summaryBy(ag_crops+Forest.coverage_bin~latitude+longitude,
                              data=full.2003.2018, FUN=c("sum"), na.rm=T)
collapsed_data <- left_join(collapsed_data.1, collapsed_data.2)
names(collapsed_data) <- c("lat","lon","zoon_outbreak_rate","avg_NDVI","fires","avg_ag_crops","avg_forest")
collapsed_data$avg_ag_crops <- ifelse(collapsed_data$avg_ag_crops>0,1,0)
collapsed_data$avg_forest <- ifelse(collapsed_data$avg_forest>0,1,0)
collapsed_data.1 <- NULL
collapsed_data.2 <- NULL


##Define countries to be plotted
african.countries <- c(
  "Algeria","Angola","Benin","Botswana","Burkina Faso","Burundi","Cabo Verde",
  "Cameroon","Central African Republic","Chad","Comoros","Democratic Republic of the Congo",
  "Republic of Congo","Ivory Coast","Djibouti","Egypt","Equatorial Guinea",
  "Eritrea","Eswatini","Ethiopia","Gabon","Gambia","Ghana","Guinea","Guinea-Bissau",
  "Kenya","Lesotho","Liberia","Libya","Madagascar","Malawi","Mali","Mauritania","Mauritius",
  "Morocco","Mozambique","Namibia","Niger","Nigeria","Rwanda","Sao Tome and Principe",
  "Senegal","Seychelles","Sierra Leone","Somalia","South Africa","South Sudan",
  "Sudan","Tanzania","Togo","Tunisia","Uganda","Zambia","Zimbabwe", "Swaziland", "Western Sahara"
)



########################################
###Plot outbreaks, ag, and forest map###
########################################

#keeping relevant data
disease <- subset(collapsed_data, zoon_outbreak_rate>0)
fordat <- subset(collapsed_data, avg_forest>0)
agzone  <- subset(collapsed_data, avg_ag_crops>0)

#Begin with a map of the African states
world <- ggplot() +
  borders("world", regions=african.countries, colour = "gray85", fill = "gray80") +
  theme_map() 

# Create the main map
map.all.col <- world +
  #ag zones
  geom_point(aes(x = lon, y = lat, alpha = avg_ag_crops),
             data = agzone, pch = 15, size = 0.8, colour = 'blue', show.legend = FALSE) +
  # forest
  geom_point(aes(x = lon, y = lat, alpha = avg_forest),
             data = fordat, pch = 15, size = 0.8, colour = 'darkgreen', show.legend = FALSE) +
  
  # Outbreaks with constant dark red color, but varying size based on outbreak rate
  geom_point(aes(x = lon, y = lat, size = zoon_outbreak_rate),
             data = disease, color = "#8B0000", alpha = .8) +  # Dark red color
  
  scale_size_continuous(name = "Average \nmonthly \noutbreak \nrates", range = c(1, 7)) +  # Adjust the size range
  
  # Move legend to the right
  theme(legend.position = c(0, 0.2),  # 0 is the far left, 0.5 centers it vertically
        legend.title = element_text(size = 10), 
        legend.text = element_text(size = 8))  # Adjust legend text size if necessary

# Save the plot as a jpeg image
jpeg("map1.jpeg", width = 6, height = 6, units = 'in', res = 500)
print(map.all.col)  # Make sure the plot is printed inside the jpeg call
dev.off()

################################
### Plot NDVI and fire map   ###
################################

# Keeping relevant data
ndvidat <- subset(collapsed_data, avg_NDVI > 0)
firedat <- subset(collapsed_data, fires >= 1)

# Remove missing coordinates
ndvidat <- na.omit(ndvidat[, c("lon", "lat", "avg_NDVI")])
firedat <- na.omit(firedat[, c("lon", "lat", "fires")])

# Begin with a map of the African states
world <- ggplot() +
  borders("world", regions = african.countries, colour = "gray85", fill = "gray80") +
  theme_map()

# Determine an appropriate tile width/height based on your data
tile_width <- 0.55  # Adjust as needed
tile_height <- 0.55  # Adjust as needed

# Create the main map
map.all.col2 <- world +
  # NDVI as background (full range) with wider tiles
  geom_tile(aes(x = lon, y = lat, fill = avg_NDVI), 
            data = ndvidat, alpha = 0.8, width = tile_width, height = tile_height) +  
  scale_fill_gradient(name = "Average NDVI", low = "white", high = "darkgreen", 
                      limits = c(0, 0.876)) +  
  
  # Fires as small circles
  geom_point(aes(x = lon, y = lat, color = fires),
             data = firedat, alpha = .8, size = 0.75) +  
  scale_color_gradient(name = "Fire Rate", low = "yellow", high = "orange", 
                       limits = c(1, max(firedat$fires, na.rm = TRUE))) +  
  
  # Adjust legend position
  theme(legend.position = "bottom",  
        legend.title = element_text(size = 10), 
        legend.text = element_text(size = 8))  

# Save the plot as a jpeg image
jpeg("map2.jpeg", width = 6, height = 6, units = 'in', res = 500)
print(map.all.col2)  # Ensure the plot is printed inside the jpeg call
dev.off()



################################
### Plot NDVI only map       ###
################################

# Keeping relevant data
ndvidat <- subset(collapsed_data, avg_NDVI > 0)

# Remove missing coordinates
ndvidat <- na.omit(ndvidat[, c("lon", "lat", "avg_NDVI")])

# Begin with a map of the African states
world <- ggplot() +
  borders("world", regions = african.countries, colour = "gray85", fill = "gray80") +
  theme_map()

# Determine an appropriate tile width/height based on your data
tile_width <- 0.55  # Adjust as needed
tile_height <- 0.55  # Adjust as needed

# Create the main map
map.all.col2a <- world +
  # NDVI as background (full range) with wider tiles
  geom_tile(aes(x = lon, y = lat, fill = avg_NDVI), 
            data = ndvidat, alpha = 0.8, width = tile_width, height = tile_height) +  
  scale_fill_gradient(name = "Average NDVI", low = "white", high = "darkgreen", 
                      limits = c(0, 0.876)) +  
  
  
  # Adjust legend position
  theme(legend.position = "bottom",  
        legend.title = element_text(size = 10), 
        legend.text = element_text(size = 8))  

# Save the plot as a jpeg image
jpeg("map3.jpeg", width = 6, height = 6, units = 'in', res = 500)
print(map.all.col2a)  # Ensure the plot is printed inside the jpeg call
dev.off()



